function [q,fq]=gp_dist(q0,f,q1)
% q0 = [1;2;3;4;5;6;7;8;9];
% f = 0.1*ones(9,1);%+[0;0;0;0.2;0;0;0;0;0];
% q1 = [1;3;5;7;9];

% q = [2;4;6;8];
% fq = [0.2;0.4;0.2;0.2];
% plot(q,f)

dq1 = q1(2) - q1(1);
q = q1(1:end-1) + dq1/2;
J = length(q1)-1;
fq = zeros(J,1);

for ii = 1:J
    fq(ii) = sum(f(q0>q1(ii) & q0<=q1(ii+1)));
end
end
